Dynamics of a deformable self-propelled domain 
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We investigate the dynamical coupling between the motion and the deformation of a single self- 
propelled domain based on two different model systems in two dimensions. One is represented 
by the set of ordinary differential equations for the center of gravity and two tensor variables 
characterizing deformations. The other is an active cell model which has an internal mechanism 
of motility and is represented by the partial differential equation for deformations. Numerical 
simulations show a rich variety of dynamics, some of which are common to the two model systems. 
The origin of the similarity and the difference is discussed. 



Introduction 

Dynamics of self-propelled objects have attracted 
much attention recently as a fundamental subject of Sta- 
tistical Physics far from equilibrium. Historically, self- 
propulsion of flexible body has been formulated in terms 
of hydrodynamics at low Reynolds number In that 
circumstance, viscous force is so large that the particle 
needs to keep non-reciprocal deformation of its shape for 
the persistent centroid motion j|2|. Swimming microor- 
ganisms are the typical examples [3, 0, S @| • Shape de- 
formation causes the center-of-mass motion in this case. 

Besides the development along this line, another class 
of self-propelled particles or domains is known in which 
the persistent motion can be maintained due to a bro- 
ken symmetry of their interfaces. In this case, interfacial 
forces are playing important roles. Experimental exam- 
ples are seen in self-propelled oil droplets in water which 
contains surfactant molecules 0, and self-propelled 
motions of vesicles in which chemical reactions take place 
0. Synthetic self-propelled systems also make a conver- 
sion of chemical energy into directed motion [HI, 
In these systems, one notes that shape deformation or 
asymmetry of chemical components around the domain is 
associated with the self-propelled motion. For example, 
an oily droplet in surfactant solution, which is spherical 
in a motionless situation becomes a banana shape when 
it undergoes a straight motion. Eukaryotic cells such as 
amoebas or fibloblast change their shape during migra- 
tion. Therefore, the coupling between the motion and the 
shape deformation is one of the most important proper- 
ties to understand the dynamics of self-propulsion from 
a unified point of view [ij, . 

From the above consideration, one may divide the self- 
propelled dynamics into two classes. One is the case that 
deformation is induced by the migration and the other 
is that the motion of the center of gravity is induced by 
the shape deformations. The oily droplets are a typical 
example of the former whereas almost all of the living 
cells belong to the latter. 

In this Letter, we consider two model systems for self- 
propelled dynamics in two dimensions. One is called a 
tensor model in terms of the velocity of the center of 
gravity and two tensor variables for deformation. The 



coupled set of equations are given by symmetry consid- 
eration and therefore they are quite general independent 
of any specific details of the self-propelled objects |15l] . 
We shall show that this model is applied, by changing 
the parameters, to both the deformation-induced motion 
and the motion-induced deformation. The other model is 
represented in the form of a partial differential equation 
for a Eucledian invariant variable of a closed loop. The 
condition of a self-propulsion is added, which is expressed 
in terms of a local deformation. Therefore this model is 
inherently a model for deformation-induced motion. By 
solving these two different model equations numerically, 
we explore possible universal and/or non- universal be- 
haviors of self-propelled dynamics. 



Deformed domain 

In this section, we introduce two model systems for a 
deformable self-propelled domain in two dimensions. One 
is based on the phenomenon of propagation of an excited 
domain in certain reaction-diffusion systems [l^. Weak 
deformation around a circular shape with radius Rq can 
be written as 



where 



R{e) = Ro + 5R{e,t) 



(1) 



(2) 



Note that since the translational motion of the domain 
will be incorporated in the velocity of the center of grav- 
ity V, the modes n — ±1 should be removed from the 
expansion ©. 

The modes z±2 represents an elliptical shape of the 
domain. We introduce a second rank tensor as <S'ii = 
—S22 = 22 + Z-2 and S12 = S21 = i{z2 - 2-2)- Similarly 
we introduce the third rank tensor from the modes z±3 
as [17| Um = 2:3 + ^-3 and U222 = ~iiz3 ~ -2-3) and 
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The time-evolution equations of v, S and U are derived 
by considering the possible couplings. Up to the third 
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order of these variables, we obtain 
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If the terms with the coefficients 02 ~ 05, &2 ~ ^5 and 
da ~ ds are ignored, these have been derived recently 
starting from the excitable reaction diffusion equations 
Furthermore, if the tensor variable U is omitted, the 
set of equations ^ and (jlj was studied previously p^. 

We shall call the dynamics of eqs. (O, Q and ^ the 
tensor model. Propagation of a domain occurs from the 
first and the second terms in eq. ([3]) for 7 > 0. The 
domain is deformed as the velocity is increased because 
of the couplings between Vi and S and U even when K2 
and K3 are positive. This case is a motion-induced de- 
formation. It should be noted that when K2 and K3 are 
negative, the domain is deformed and causes a drift mo- 
tion due to the couplings SU and SUU in eq. ([3]). This 
implies a deformation-induced motion. 

Another model for a deformation-induced motion is an 
active cell model which is expressed in term of the partial 
differential equation for a closed domain boundary X(s) 
for < s < L in two dimensions where L is the boundary 
length and X(0) — X(L). Here we employ the intrinsic 
representation of a closed loop as (l9| 



dXjs) 
ds 



t(s) , 



(6) 



with t(s) the tangential unit vector. The Frenet-Serret 
formula gives us 



dt{s) 
ds 



K(s)n(s) , 



(7) 



where n is the curvature and n is the unit normal which 
is written as 

2tt 27r 
n(s) = (cos— (<?!>(s) -f-s),sin— ((/.(s) -l-s)) , (8) 

where (j){s) represents deformation around a circular 
shape. Throughout this paper, we assume that is suf- 
ficiently small and is a single valued function of s. 

As a phenomenological description of the active dy- 
namics of (j), we make a symmetry argument. Because 
of the isotropy of space and the parity symmetry, the 
dynamics of </> should be invariant against the following 
transformations: (i) + , (ii)s s + sq with 00 
and So constants and (iii)s — s and (f> —0. Keeping 
these in mind, we write down the equation for (j) up to 
bilinear order of 6 as 



34(9.0) (920) - 



(9) 



Here we assume that there is no instability in the long 
wave length limit so that gi is non-negative. To make 
the circular shape is unstable, then, we impose that §2 is 
positive and 173 is also positive to recover the stability at 
the short wave length region. Under this condition the 
term with g\ is expected to be irrelevant and can be ig- 
nored. The nonlinear term with the coefficient g^ can be 
written in a variational form 5 J ds{ds(j>)^ /Scj). (To make 
the potential functional bounded below, we need to add 
J ds{ds(f>)'^ ■) Therefore, this nonlinearity does not cause 
any asymptotic complex dynamics such as domain oscil- 
lation and should be ignored. The last term is not vari- 
ational. Under these considerations, the minimal non- 
trivial equation for is given by 



(10) 



where the coefficients are eliminated by redefining t, s 
and 0. The sign in front of the last term is chosen to be 
negative without loss of generality. Note that the param- 
eter which we can control is only the domain boundary 
length L. In a previous paper, eq. (jlOp was derived 
approximately starting from the free energy functional 
for the interfacial energy and the curvature energy of a 
domain [lH. 

Equation (jlOp describes the dynamics of deformations. 
By using eqs. ^ and (O, the shape of the do- 

main is determined. It should be noted, however, that 
the translational motion cannot be obtained by the so- 
lution. We have to impose the condition for the time- 
dependence of the center of mass. Since we are consider- 
ing a deformation-induced motion, it should depend on 
the curvature. By taking account of the fact that the 
normal unit n is the basic vector variable, the velocity of 
the center of gravity should take the following form 
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where Y is an unknown function of k. Since the defor- 
mation is weak as we have assumed, we may expand Y 
in powers of k as F = ao + aiK + a2K^ + a3{dsK)^ ■■■■ The 
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lowest order term vanishes because 

rL 

dsn = . 



(12) 







This is the condition that the boundary is closed. The 
first order term also vanishes identically because of eq. 

0; 



dsnn — 



As a result, the velocity is given by 



(13) 
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and 



where we have used the relations k = (2n/ L){l+ds 
eqs. (fT2|) and (fT3|) . Note that the terms with the coeffi- 
cients a2 and in eq. (jl4p correspond to the terms with 
(74 and 175 of eq. ([9]) respectively. Since we have ignored 
the 54 term, we retain only the as term for consistency 



(2^)^a3 
L3 



ds(a20)2n 



(15) 



Equations pUj) and ([T5|) complete the motion of a do- 
main. 



Numerical simulations I 

First we show the results of numerical simulations of 
the reduced tensor model for the coupled set of equa- 
tions for Vi and Sij ignoring Uijk- This is justified when 
the relaxation of U is sufficiently rapid, i.e., K3 is large 
enough and the velocity is sufiiciently small. That is, we 
consider eqs. ([3]) and (|3]) with 02 = 03 = 05 = ae = 62 = 
^5 = ^6 = 0. The terms with the coefficients 04 and 64 
are also omitted. Furthermore, we allow the case that 
K2 is negative in eq. In this section, we put 63 = 1 
without loss of generality. 

The set of equations has been solved numerically for 
ai = 1 and bi = 0.5 and changing the parameters K2 and 
7. The simple Euler scheme has been employed with the 
time increment At = 10"^. We have checked the numer- 
ical accuracy by using At — 10^*. The phase diagram 
obtained is displayed in Fig. [1] In the region indicated 
by the cross symbol, the domain is motionless whereas 
it undergoes a straight motion in the region of the open 
squares and a circular motion in the region of the open 
circles. These are essentially the same as the previous 
findings [l^. The solid line is the boundary between 
the motionless state and the straight motion whereas the 
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FIG. 1 Phase diagram obtained numerically for eqs. ([S]) and 
(|4| ignoring the Uijk terms. The parameters are set to be 
a\ = 1.0 and &i = 0.5. The meanings of the synbols and the 
lines are given in the text. 



broken line is the phase boundary between the straight 
motion and the circular motion. The dotted line indicates 
the subcritical hopf-bifurcation line from the straight mo- 
tion to the rectangular motion which is described below. 
This line has been obtained as the stability limit of the 
straight motion. 

Two interesting dynamics appear for K2 < 0. One is 
the so called rectangular motion in which the domain 
repeats a straight motion and stopping alternatively as 
shown in Fig. [5] for 7 = —0.04 and K2 = —0.1. This oc- 
curs in the region indicated by the solid squares in Fig. 
[1] During the stopping interval the domain changes the 
shape and the propagation direction almost by 90°. Ei- 
ther clock-wise rotation or counter clock-wise rotation 
seem to occur at random and may depend on noises 
caused unavoidably in the numerical computations. In 
the most of the region shown by the solid circles, where 
7 > and K2 < 0, a kind of circular motion is observed. 
However, this circular motion does not have a single fre- 
quency but has a multi-frequency with irrational ratio 
and therefore the motion is quasi-periodic as shown in 
Fig. [3Ja) for 7 = 0.02 and K2 = —0.06. In order to 
confirm that the motion is quasi-periodic, we have ana- 
lyzed the return map as shown in Fig. [Sfb) where the 
values of the x-componcnt of the location of the domain 
are plotted every time that the domain crosses the line 
—3 < a; < 3 and y — 3. 



Numerical simulations II 

The full set of equations (O, (g]) and (P (but with 
the simplification 02 = 03 = 04 = 05 = ag = 64 = 
— be = ds = di = d^ — de — have also been 
solved numerically. The modified Euler method with the 
time increment either At = 10"^ or At = 10~* has been 
employed. The parameters are chosen as ai = —1.0, 
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FIG. 3 (a) Quasi-periodic motion for a — 1.0, b = 0.5, H2 = 
FIG. 2 Rectangular motion for ai = 1.0, bi = 0.5, K2 = -0.06 and 7 = 0.02 rotating in the counterclockwise direction. 
-0.1 and 7 = -0.04. The arrows and the digits indicate (b) Return map of the x coordinate, 
the direction of motion and the time-sequence of the motion 
respectively. 



61 = -0.5, 62 = 0.3, di = 0.1, d2 = 0.8 and 7 = 1.0. 
The relaxation rates K2 and K3 are varied. Since these 
relaxation rates are chosen to be positive in this section, 
the cubic term in eqs. (jl]) and dU are not considered, 
i.e., 63 = = 0. 

The phase diagram is obtained as shown in Fig. Ufa). 
The straight motion and the circular motion appear in 
the region indicated by the squares and by the circles re- 
spectively. In the region indicated by the stars for the 
smaller values of K2 and for K3 = 0.1, the domain motion 
becomes chaotic. In order to confirm the chaotic behav- 
ior, we have evaluated the maximum Lyapunov exponent, 
Ai, associate with the domain trajectory as depicted in 
Fig. m^b). It is evident that Ai becomes positive for 
K2 < 0.8. For smaller values of K2 we encounter a nu- 
merical instability and cannot obtain any accurate value 
of the exponent. 

When K2 is large and K3 is small, a zig-zag motion is ob- 
served as indicated by the triangles. The time-sequence 
of the snapshots for K2 = 0.9 and K3 =0.1 is displayed 
in Fig. [Sla). The domain is traveling from the left to 
the right. The angle of the zig-zag motion is about 60°. 
A chaotic trajectory is displayed in Fig. ^h) which is 
obtained for K2 = 0.5 and K3 = 0.1. 



Numerical simulations III 

The tensor model takes account only of two long wave- 
length deformation modes. In contrast, there is no re- 
striction in the active cell model [l^ because the dy- 
namics is represented by the partial differential equation 
pn)) . We use the spectral method with the 4th order 
Runge-Kutta algorithm with the total number of wave 
modes N = 128 to solve eq. dTU]) together with ([T5|). 
The geometrical condition must be satisfied at each 
time step, which is represented in terms of the complex 
variable ni + in2 = exp[(27ri/_L)(s -I- (^(s))] as 




ds exp[{2Tri/L){s + (/)(s))] = . (16) 



However, since the model ignores the Oi(t>^) terms, this 
condition is not automatically fulfilled [1^, [l^ . In order 
to overcome this problem, we employ the Bayesian esti- 
mate First we replace the integrand of eq. by 
exp[(27ri/L)(s -|- (j){s) + e(s))] where e(s) is the unknown 
error function. We introduce the likelihood function as 
well as the a priori distribution for the modes of e(s). By 
optimizing the logarithmic likelihood, we determine e(s) 
to enclose the domain boundary. 

For an internal consistency, we have verified numeri- 
cally that the length L is actually unchanged appreciably 
by these numerical methods. It is noted that the domain 
area is time-dependent in this model system. 

As is seen in Figs. [T]and|4l the tensor model has exhib- 
ited several types of deformation dynamics by changing 
the three parameters 7, K2 and K3. In contrast, the active 
cell model produces these similar dynamics changing only 
the value of L. Numerical simulations of Eqs. (ITOl) and 
psp show that the intricate dynamics appear in the three 
characteristic windows; Wi ~ {L\L e [12.0, 14.4]}, W2 ~ 
{L\L e [21.0,21.9]}, and W3 ~ {L\L G [29.0,29.9]}. Cir- 
cula r, q uasi-periodic, and rectangular motions emerge in 
Wi [18'], quasi-periodic and chaotic motions in W2, and 
zig-zag motions in W3 . These motions occur successively 
by changing the value of L. The straight motion is ob- 
tained in the region 6.2 < L < 12.0. No shape insta- 
bility occurs for L < 6.2 and therefore no propagation 
of domain. There is an obvious reason why the complex 
dynamics appear in the windows Wi, W2 and W3. The 
linear stability analysis of Eqs. pO)) about the trivial 
solution ^ = shows the Fourier mode-n deformation 
becomes unstable at L = 2?™. Actually the mode-2 be- 
comes unstable in Wi = {L\L ^ 47r}, the mode-3 in 
W2 = {L\L - 6n} and the mode-4 in W3 = {L\L - 87r}. 
Since the codimension two bifurcation points exist near 
these critical points, the various types of motion would 
appear. 

Figure M,^) displays the trajectory obtained in W2, 
where an apparently chaotic motion appears. The broad 
power spectrum of the Fourier amplitudes of (/> shown 
in the inset implies that the time-evolution of (j){s,t) is 
a kind of spatio-temporal chaos. Figure \^h) shows a 
quasi-periodic solution obtained in W2. Figures [Hllc) and 
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FIG. 4 (a) Phase diagram for the full set of equations for ai — 
-1.0, bi = -0.5, b2 = 0.3, di = 0.1, d2 = 0.8 and 7 = 1.0. 
The meaning of the symbols is given in the text. In the region 
indicated by the symbol x , an numerical instability occurs so 
that we have no definite conclusion about the motion, (b) 
Maximum Lyapunov exponent obtained numerically for K3 = 
0.1. The other parameters are the same as those of the phase 
diagram (a). 



(d) show the two types of trajectory obtained in W3. Al- 
though the trajectory in Fig. [SJc) seems complicated, 
the orderly turns occurring in the trajectory (@, (3), in 
Fig. M.c)) have almost 180° ± 45°. Thus this is a kind 
of zig-zag motion, however unlike the zig-zag motion in 
Fig. EI a), the effective modes governing the deformation 
are n = ±4 here. The trajectory in Fig. [DJd) is also con- 
sidered to be a kind of zig-zag motion, because the main 
part of the trajectory indicated by the dotted gray line 
shows ±45° turns. The zig-zag motion with n — ±3 as 
in Fig, [Sfa) has not been found in the active cell model. 

Coupled set of equations for the Fourier amplitudes 

The set of equations dU and ([5]) can be writ- 
ten in terms of the Fourier components in eq. ©. 
Here we define the complex variables as zi = vi — iv2, 
Z2 = ^{Sii - iS'12) and Z3 = ^(C/in + ^^7222)• The time- 
evolution equations for zi, Z2 and Z3 are given from eqs. 
®, Q and ® by 

zi = (7 + 0^111211^ + ^12122^ + ^131231^)21 

-t- eiiiiZ2 + ei22"i^23 + 6132223 + 6142322 , (17) 

22 = (-K2 + rf2l|2ip + rf22|22p + (i23|23|^)22 

+ 62121 -I- e222"l23 -|- e232"2232l , (18) 

23 = (-K3 + rf3l|2lp -t- (i32|22p H- d33|23p)23 

-I- esizf + 6322122 + 6332"l2| , (19) 

where the dot means the time derivative and the bar 
indicates the complex conjugate. All the coefficients are 
real and are given by dn = —1, rfi2 = —804, rfi3 = 

-1605, d21 = -h, d22 = -863, ^23 = -I665, ^31 = -d4, 

d32 = -Sds, d33 = -16(^3, 611 = -2ai, 612 = -2a2, 
613 = -8a3, 614 IGflg, 621 = 5i/4, 622 = b2, 623 = 266, 
631 = di/8, 632 = c?2/2 and 633 — 2dg. It is important 
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FIG. 5 (a) Zigzag motion for K2 = 0.9 and K3 = 0.1. The 
domain moves from the left to the right, (b) Chaotic motion 
for K2 = 0.5 and K3 — 0.1. 



to note that this set of equations is invariant under the 
transformation 

(21, 22, 23) ^ {e''zi,e^'''z2,e^'''z3) , (20) 

for an arbitrary phase angle 9. This invariance arises 
from the isotropy of space. 

When the variable 23 is omitted, the set of equations 
(HH) and (UHl) is the same as those considered by Arm- 
bruster et al They motivated to study some par- 

tial differential equation like the Kuramoto-Sivashinsky 
equation in one dimension under the periodic bound- 
ary condition [2l| and had no consideration of the self- 
propelled domain dynamics. In fact, there are the follow- 
ing correspondences; the motionless circular shape do- 
main o the trivial solution, the deformed motionless do- 
main o pure mode, the straight motion o the standing 
wave, the rectangular motion i-)- the heteroclinic cycle, 
the rotating motion o the traveling wave and the quasi- 
periodic motion the modulated wave. The former is 
the motions obtained in our theory whereas the latter is 
the terminology of Armbruster et al [20| . 

The variable cj) in the active cell model can also be 
represented in terms of the Fourier modes. If the modes 
higher than n — S relax rapidly to the stationary values, 
we may retain only the modes n = 1,2,3. The time- 
evolution equations are essentially the same structures 
as those in eqs. ([T7]), and P^ . Consequences of 
this will be discussed below. 



Discussion 

We have shown that various self-propelled motions 
appear both in the tensor model and the active cell 
model. The dynamics common to these two systems are 
straight motion, circular periodic and quasi-periodic mo- 
tions, rectangular motion, zig-zag motions and chaotic 
motion. 

Here we note the main difference between our 
model and the self-propelled swimmers at low Reynolds 
number[2|- Our model equations are autonomous, thus 
the shape deformation and the centroid migration are 
spontaneously created. On the other hand in the latter 
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FIG. 6 Trajectories of a domain in the active cell model for 
as = -4 X W^L. (a) Chaotic motion for L = 21.6. The 
inset shows the power spectrum calculated from the defor- 
mation modes, (b) Quasi-periodic motion for L = 21.8. (c) 
Zig-zag motion for L — 29.2 with the turn angle 180° ± 45°. 
The circled numbers mean the qualifying order of trajectory, 
(d) Zig-zag motions for L = 29.6. The "coarse-grained" tra- 
jectory (the dotted gray line) shows ±45° turns. 



frameworks, the deformation of flexible body is opera- 
tionally given and resulting motion is considered within 
Stokes dynamics. Owing to the autonomous proper- 
ties, our models exhibit successive bifurcations leading 
to richer dynamics. Since the tensor model ([3]), (|H), and 
([5]) have been derived from the excitable reaction diffu- 
sion model, the bifurcation from a simple to a complex 
motion should be explored experimentally in physico- 
chemical systems such as oily droplet systems. Along this 
line, propagating actin waves and recovery of actin poly- 
merization from complete depolymerization observed in 
Dictyostelium cells [23| might give a clue to a connection 
between our model and possible bifurcations in the dy- 
namics of living cells. 

Moreover, we note the difference of the zigzag motions 
in the two models. The zig-zag motion with the angle 
about 60° has been obtained in the tensor model. This is 
attributed to the fact that only the second and the third 
modes are considered in the tensor model so that the 
deformation with three fold symmetry is possible which 
triggers the 60° zigzag motion. In the active cell model, 
on the other hand, more complicated zigzag motion with 
the angle about 135° and 45° appears as in Figs. [6jc) and 
(d) where L/2tt is close to 4. Therefore, it is expected 



that higher modes such as the fourth mode are dominant 
for the zig-zag motion in the active cell model. 

We emphasize that the rectangular motion, the 60° 
zigzag motion and an apparently chaotic motion have 
been observed in real experiments of amoebas ^J, IT3 |. 
Therefore the present approach based on the symmetry 
argument to construct the time-evolution equations cap- 
tures the essential feature of the coupling between the 
shape and the motion of a self-propelled domain. 
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